% Test the function makeppft filter.
%
% Yoel Shkolnisky, December 2010.
global verbose
for precond=0:1
    for n=2:2:128
        A=rand(n,n);

        tic;
        [pp1 pp2]=optimizedPPFT2D(A);
        if precond
            B=preconditionedAdjointPPFT2D(pp1,pp2);
        else
            B=optimizedAdjointPPFT2D(pp1,pp2);
        end
        t1 = toc;
        L=3*n-2; %n+2n-1-1

        filename=sprintf('filters\\Hhatfilter%dprecond%d.dat',n,precond);
        fid=fopen(filename,'r');
        if fid>0
            N=fread(fid,1,'int');
            timestamp=fread(fid,1,'real*8');
            if verbose
                fprintf('loading 2D Hhat filter which was created on %s for n = %d and precond %d\n',datestr(timestamp),n,precond);
            end
            Hhat=fread(fid,round(L*(L+1)/2),'real*8');
            Hhat=Hhat+1i*fread(fid,round(L*(L+1)/2),'real*8');
            Hhat=reshape(Hhat,1,round(L*(L+1)/2));
            fclose(fid);
        else
            timestamp=now;
            if verbose
                fprintf('creating 2D Hhat filter on %s with n = %d and precond %d for future usage\n',datestr(timestamp),n,precond);
            end
            H=makeHFilter2D(n,precond);
            Hhat=fftn(H.filter,[L L]);
            fid=fopen(filename,'w');
            if fid<0
                error('failed to open file for writing, filter was not created\n');
            else
                fwrite(fid,n,'int');
                fwrite(fid,timestamp,'real*8');
                for i = 1:L
                    fwrite(fid,real(Hhat(i,i:L)),'real*8');
                end
                for i = 1:L
                    fwrite(fid,imag(Hhat(i,i:L)),'real*8');
                end
                fclose(fid);
                fid=fopen(filename,'r');
                if fid>0
                    N=fread(fid,1,'int');
                    timestamp=fread(fid,1,'real*8');
                    if verbose
                        fprintf('loading 2D Hhat filter which was created on %s for n = %d and precond %d\n',datestr(timestamp),n,precond);
                    end
                    Hhat=fread(fid,round(L*(L+1)/2),'real*8');
                    Hhat=Hhat+1i*fread(fid,round(L*(L+1)/2),'real*8');
                    Hhat=reshape(Hhat,1,round(L*(L+1)/2));
                    fclose(fid);
                end
            end
        end

        idx=[];
        jdx=1:L;
        kdx=L:-1:1;

        tic;
        Ahat=fftn(A,[L L]);

        for i = 1:L
            Htmp=[Hhat(idx) Hhat(jdx)];
            if i > 1
                idx=[i+1 idx+kdx(1:i-1)];
            else
                idx=i+1;
            end
            jdx=jdx(2:end)+L-i;
            Ahat(i,:)=Ahat(i,:).*Htmp;
        end
        Cc=ifftn(Ahat,[L L]);
        Cc=Cc(n:2*n-1,n:2*n-1);
        t2=toc;

        reportTestResult(n,(norm(B(:)-Cc(:))./norm(B(:))),t1/t2);    
    end
end
